
NUM_MC = 10000
set.seed(0)

# SIMULATION PARAMETERS ################################################
N = 75
T = 75
b = matrix(c(0.1,0.2,0.1,0.1,0.1),,1)
DGP_AR = TRUE ### AR if TRUE; MA if FALSE ###
rho = .25 ### In case of AR ###
q = 10
MAcoef = matrix(c(rep(1/q,q)),,1) ### In case of MA ###
true_beta = matrix(c(.1,.1),2,1)

# TUNING FOR THOMPSON ##################################################
m_thompson = 2

# MENZEL BOOTSTRAP ITERATIONS ##########################################
num_boot = 1000
# MNW BOOTSTRAP ITERATIONS #############################################
num_boot_MNW = 1000

# MAMEN RANDOM VARIABLE ################################################
mammen <- function(para){ a = para[1]; b = para[2]; p = para[3]; return( (p*a+(1-p)*b)^2 + (p*a^2+(1-p)*b^2-1)^2 + (p*a^3+(1-p)*b^3-1)^2 ); }
mammen_para = nlm(mammen,c(-1,1,0.5))$estimate
rmammen <- function(n,para){(para[1]-para[2])*(runif(n)<para[3])+para[2]}

# FUNCTION: EIGENVALUE-CORRECTION ######################################
EV_Correction <- function( A ){
 EV = eigen(A)
 Eval = EV$values
 Eval = (Eval >= 0) * Eval + (Eval < 0 ) * 0
 Evec = EV$vectors
 return( (Evec) %*% diag(Eval) %*% t(Evec) )
}

########################################################################
# BEGIN MONTE CARLO ####################################################
########################################################################
list_IID_95 = NULL
list_PNW_95 = NULL
list_CRi_95 = NULL
list_CRt_95 = NULL
list_CGM_95 = NULL
list_MNW_95 = NULL
list_MEN_95 = NULL
list_ME2_95 = NULL
list_THO_95 = NULL
list_DDG_95 = NULL
list_CHS_95 = NULL
for( iter in 1:NUM_MC ){

########################################################################
# DATA GENERATION (BEGIN)
########################################################################

### ALPHA (NORMAL) (N by 1) ############################################
alpha_0 = matrix(rnorm(N,0,1),N,1)
alpha_1 = matrix(rnorm(N,0,1),N,1)
alpha_2 = matrix(rnorm(N,0,1),N,1)
alpha_3 = matrix(rnorm(N,0,1),N,1)

### GAMMA (AR NORMAL) (1 by T) #########################################
if( DGP_AR ){
 gamma_0 = matrix(rnorm(1,0,1),1,T)
 gamma_1 = matrix(rnorm(1,0,1),1,T)
 gamma_2 = matrix(rnorm(1,0,1),1,T)
 gamma_3 = matrix(rnorm(1,0,1),1,T)
 for( idx in 2:T ){
  gamma_0[1,idx] = rho*gamma_0[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_1[1,idx] = rho*gamma_1[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_2[1,idx] = rho*gamma_2[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_3[1,idx] = rho*gamma_3[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
 }
}

### GAMMA (MA NORMAL) (1 by T) #########################################
if( !DGP_AR ){
 pre_gamma_0 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_1 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_2 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_3 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 gamma_0 = matrix(0,1,T)
 gamma_1 = matrix(0,1,T)
 gamma_2 = matrix(0,1,T)
 gamma_3 = matrix(0,1,T)
 for( idx in 1:T ){
  gamma_0[1,idx] = t(MAcoef) %*% matrix(pre_gamma_0[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_1[1,idx] = t(MAcoef) %*% matrix(pre_gamma_1[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_2[1,idx] = t(MAcoef) %*% matrix(pre_gamma_2[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_3[1,idx] = t(MAcoef) %*% matrix(pre_gamma_3[1:(dim(MAcoef)[1])+idx-1],,1)
 }
}

### EPSILON (NORMAL) (N by T) ##########################################
epsilon_0 = matrix(rnorm(N*T,0,1),N,T)
epsilon_1 = matrix(rnorm(N*T,0,1),N,T)

X = b[3,1] * ( alpha_1 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_2 ) +
    b[4,1] * ( alpha_2 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_1 ) +
    b[5,1] * epsilon_0

U = b[1,1] * ( alpha_0 %*% matrix(1,1,T) ) +
    b[2,1] * ( matrix(1,N,1) %*% gamma_0 ) +
    b[3,1] * ( alpha_1 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_3 ) +
    b[4,1] * ( alpha_3 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_1 ) +
    b[5,1] * epsilon_1

Y = true_beta[1,1] + true_beta[2,1] * X + U

### DOUBLE DIFFERENCE Y, X, AND U ######################################
Y = Y - matrix( rowMeans(Y), N, T ) - t( matrix( colMeans(Y), T, N ) ) + matrix( mean(Y), N, T )
X = X - matrix( rowMeans(X), N, T ) - t( matrix( colMeans(X), T, N ) ) + matrix( mean(X), N, T )
U = U - matrix( rowMeans(U), N, T ) - t( matrix( colMeans(U), T, N ) ) + matrix( mean(U), N, T )

########################################################################
# DATA GENERATION (ENDED)
########################################################################



########################################################################
# ESTIMATION (BEGIN)
########################################################################

### OLS ################################################################
beta_hat = solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% t(cbind(1,matrix(c(X),,1))) %*% matrix(c(Y),,1)

### RESIDUALS ##########################################################
U_hat = Y - beta_hat[1,1] - beta_hat[2,1] * X

########################################################################
# ESTIMATION (ENDED)
########################################################################



########################################################################
# ESTIMATE IID OMEGA (BEGIN)
########################################################################
Omega_hat_iid = matrix(0,2,2)
Omega_hat_iid[1,1] = 0 # t(matrix(c(1*U_hat),,1)) %*% matrix(c(1*U_hat),,1)
Omega_hat_iid[2,1] = 0 # t(matrix(c(X*U_hat),,1)) %*% matrix(c(1*U_hat),,1)
Omega_hat_iid[1,2] = 0 # t(matrix(c(1*U_hat),,1)) %*% matrix(c(X*U_hat),,1)
Omega_hat_iid[2,2] = t(matrix(c(X*U_hat),,1)) %*% matrix(c(X*U_hat),,1)
Omega_hat_iid = 1/(N*T) * Omega_hat_iid
########################################################################
# ESTIMATE IID OMEGA (ENDED)
########################################################################



########################################################################
# MENZEL'S BOOTSTRAP (BEGIN)
########################################################################

Ubar = mean(U_hat)
XUbar = mean(X*U_hat)
a1 = rowMeans(U_hat)
a2 = rowMeans(X*U_hat)
g1 = colMeans(U_hat)
g2 = colMeans(X*U_hat)
w1 = U_hat - matrix(a1,N,T) - t(matrix(g1,T,N))
w2 = X*U_hat - matrix(a2,N,T) - t(matrix(g2,T,N))

S2_a1 = rowSums((U_hat - Ubar)^2)/(N-1)
S2_a2 = rowSums((X*U_hat - XUbar)^2)/(N-1)
S2_g1 = colSums((U_hat - Ubar)^2)/(T-1)
S2_g2 = colSums((X*U_hat - XUbar)^2)/(T-1)
S2_w1 = sum((U_hat- matrix(a1,N,T) - t(matrix(g1,T,N)) - Ubar)^2) / (N*T - N - T)
S2_w2 = sum((X*U_hat - matrix(a2,N,T) - t(matrix(g2,T,N)) - XUbar)^2) / (N*T - N - T)

sigma2_a1 = max(c(0,S2_a1 - S2_w1/T))
sigma2_a2 = max(c(0,S2_a2 - S2_w2/T))
sigma2_g1 = max(c(0,S2_g1 - S2_w1/N))
sigma2_g2 = max(c(0,S2_g2 - S2_w2/N))
sigma2_w1 = S2_w1
sigma2_w2 = S2_w2

eps_hat = U_hat - matrix(a1,N,T) - t(matrix(g1,T,N)) - mean(U_hat)
mu4_e = sqrt(2)*max(c(0.1,sqrt(mean(eps_hat^4)-mean(eps_hat)^4)))
kappa_a = mu4_e * log(T)
kappa_g = mu4_e * log(N)

D_a1 = 1*(T*sigma2_a1 >= kappa_a)
D_a2 = 1*(T*sigma2_a2 >= kappa_a)
D_g1 = 1*(N*sigma2_g1 >= kappa_g)
D_g2 = 1*(N*sigma2_g2 >= kappa_g)

lambda_a1 = (D_a1 * T * sigma2_a1) / (D_a1*T*sigma2_a1 + sigma2_w1)
lambda_a2 = (D_a2 * T * sigma2_a2) / (D_a2*T*sigma2_a2 + sigma2_w2)
lambda_g1 = (D_g1 * N * sigma2_g1) / (D_g1*N*sigma2_g1 + sigma2_w1)
lambda_g2 = (D_g2 * N * sigma2_g2) / (D_g2*N*sigma2_g2 + sigma2_w2)

lambda_a1_No_Model_Sel = (1 * T * sigma2_a1) / (1*T*sigma2_a1 + sigma2_w1)
lambda_a2_No_Model_Sel = (1 * T * sigma2_a2) / (1*T*sigma2_a2 + sigma2_w2)
lambda_g1_No_Model_Sel = (1 * N * sigma2_g1) / (1*N*sigma2_g1 + sigma2_w1)
lambda_g2_No_Model_Sel = (1 * N * sigma2_g2) / (1*N*sigma2_g2 + sigma2_w2)

boot_beta = NULL
for ( boot_iter in 1:num_boot ){

# DRAW BOOTSTRAP SAMPLES ###############################################
ki = sample(1:N,N,replace=TRUE)
kt = sample(1:T,T,replace=TRUE)

# STEP (a) #############################################################
a1_star = a1[ki]
a2_star = a2[ki]
g1_star = g1[kt]
g2_star = g2[kt]
# STEP (b) #############################################################
omega1_star = matrix(rmammen(N,mammen_para),N,T)
omega2_star = t(matrix(rmammen(T,mammen_para),T,N))
w1_star = omega1_star*omega2_star*w1[ki,kt]
w2_star = omega1_star*omega2_star*w2[ki,kt]
# STEP (c) #############################################################
z1_star = lambda_a1^0.5*matrix(a1_star,N,T) + lambda_g1^0.5*t(matrix(g1_star,T,N)) + w1_star
z2_star = lambda_a2^0.5*matrix(a2_star,N,T) + lambda_g2^0.5*t(matrix(g2_star,T,N)) + w2_star
z1_star_No_Model_Sel = lambda_a1_No_Model_Sel^0.5*matrix(a1_star,N,T) + lambda_g1_No_Model_Sel^0.5*t(matrix(g1_star,T,N)) + w1_star
z2_star_No_Model_Sel = lambda_a2_No_Model_Sel^0.5*matrix(a2_star,N,T) + lambda_g2_No_Model_Sel^0.5*t(matrix(g2_star,T,N)) + w2_star
# STEP (d) #############################################################
boot_beta = cbind(boot_beta, beta_hat + N*T*solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% matrix(c(mean(z1_star),mean(z2_star)),2,1))
boot_beta_No_Model_Sel = cbind(boot_beta, beta_hat + N*T*solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% matrix(c(mean(z1_star_No_Model_Sel),mean(z2_star_No_Model_Sel)),2,1))
}

cover_MEN = abs(beta_hat-true_beta) < apply(boot_beta,1,sd)*qnorm(0.975)
list_MEN_95 = rbind(list_MEN_95, t(cover_MEN))
cover_ME2 = abs(beta_hat-true_beta) < apply(boot_beta_No_Model_Sel,1,sd)*qnorm(0.975)
list_ME2_95 = rbind(list_ME2_95, t(cover_ME2))

########################################################################
# MENZEL'S BOOTSTRAP (ENDED)
########################################################################



########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1 = matrix(0,2,2)
for( i in 1:N ){
 # Omega_hat_1[1,1] = Omega_hat_1[1,1] + sum( t(matrix(    1*U_hat[i,],1,T)) %*% matrix(    1*U_hat[i,],1,T) )
 # Omega_hat_1[2,1] = Omega_hat_1[2,1] + sum( t(matrix(X[i,]*U_hat[i,],1,T)) %*% matrix(    1*U_hat[i,],1,T) )
 # Omega_hat_1[1,2] = Omega_hat_1[1,2] + sum( t(matrix(    1*U_hat[i,],1,T)) %*% matrix(X[i,]*U_hat[i,],1,T) )
 Omega_hat_1[2,2] = Omega_hat_1[2,2] + sum( t(matrix(X[i,]*U_hat[i,],1,T)) %*% matrix(X[i,]*U_hat[i,],1,T) )
}
Omega_hat_1 = min(c(N,T))/N/(N*T^2) * Omega_hat_1

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2 = matrix(0,2,2)
for( t in 1:T ){
 # Omega_hat_2[1,1] = Omega_hat_2[1,1] + sum( matrix(    1*U_hat[,t],N,1) %*% t(matrix(    1*U_hat[,t],N,1)) )
 # Omega_hat_2[2,1] = Omega_hat_2[2,1] + sum( matrix(X[,t]*U_hat[,t],N,1) %*% t(matrix(    1*U_hat[,t],N,1)) )
 # Omega_hat_2[1,2] = Omega_hat_2[1,2] + sum( matrix(    1*U_hat[,t],N,1) %*% t(matrix(X[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[2,2] = Omega_hat_2[2,2] + sum( matrix(X[,t]*U_hat[,t],N,1) %*% t(matrix(X[,t]*U_hat[,t],N,1)) )
}
Omega_hat_2 = min(c(N,T))/T/(N^2*T) * Omega_hat_2

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3 = matrix(0,2,2)
# Omega_hat_3[1,1] = t(matrix(c(1*U_hat),,1)) %*% matrix(1*c(1*U_hat),,1)
# Omega_hat_3[2,1] = t(matrix(c(X*U_hat),,1)) %*% matrix(1*c(1*U_hat),,1)
# Omega_hat_3[1,2] = t(matrix(c(1*U_hat),,1)) %*% matrix(1*c(X*U_hat),,1)
Omega_hat_3[2,2] = t(matrix(c(X*U_hat),,1)) %*% matrix(1*c(X*U_hat),,1)
Omega_hat_3 = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3

### OMEGA HAT 4TH TERM FOR THOMPSON ####################################
Omega_hat_4_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_4_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  # temp_Omega_hat_4_thompson[1,1] = temp_Omega_hat_4_thompson[1,1] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
  # temp_Omega_hat_4_thompson[2,1] = temp_Omega_hat_4_thompson[2,1] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
  # temp_Omega_hat_4_thompson[1,2] = temp_Omega_hat_4_thompson[1,2] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,2] = temp_Omega_hat_4_thompson[2,2] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_4_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_4_thompson
 Omega_hat_4_thompson = Omega_hat_4_thompson + temp_Omega_hat_4_thompson
}
Omega_hat_4_thompson = min(c(N,T))/T * Omega_hat_4_thompson

### OMEGA HAT 5TH TERM FOR THOMPSON ####################################
Omega_hat_5_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_5_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  # temp_Omega_hat_5_thompson[1,1] = temp_Omega_hat_5_thompson[1,1] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
  # temp_Omega_hat_5_thompson[2,1] = temp_Omega_hat_5_thompson[2,1] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
  # temp_Omega_hat_5_thompson[1,2] = temp_Omega_hat_5_thompson[1,2] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,2] = temp_Omega_hat_5_thompson[2,2] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
 }
 temp_Omega_hat_5_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_5_thompson
 Omega_hat_5_thompson = Omega_hat_5_thompson + temp_Omega_hat_5_thompson
}
Omega_hat_5_thompson = min(c(N,T))/T * Omega_hat_5_thompson

### OMEGA HAT 6TH TERM FOR THOMPSON ####################################
Omega_hat_6_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_6_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  # temp_Omega_hat_6_thompson[1,1] = temp_Omega_hat_6_thompson[1,1] + sum( matrix(      1*U_hat[,t],N,1) * (matrix(      1*U_hat[,t-j],N,1)) )
  # temp_Omega_hat_6_thompson[2,1] = temp_Omega_hat_6_thompson[2,1] + sum( matrix(  X[,t]*U_hat[,t],N,1) * (matrix(      1*U_hat[,t-j],N,1)) )
  # temp_Omega_hat_6_thompson[1,2] = temp_Omega_hat_6_thompson[1,2] + sum( matrix(      1*U_hat[,t],N,1) * (matrix(X[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,2] = temp_Omega_hat_6_thompson[2,2] + sum( matrix(  X[,t]*U_hat[,t],N,1) * (matrix(X[,t-j]*U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_6_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_6_thompson
 Omega_hat_6_thompson = Omega_hat_4_thompson + temp_Omega_hat_6_thompson
}
Omega_hat_6_thompson = min(c(N,T))/T * Omega_hat_6_thompson

### OMEGA HAT 7TH TERM FOR THOMPSON ####################################
Omega_hat_7_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_7_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  # temp_Omega_hat_7_thompson[1,1] = temp_Omega_hat_7_thompson[1,1] + sum( matrix(      1*U_hat[,t-j],N,1) * (matrix(      1*U_hat[,t],N,1)) )
  # temp_Omega_hat_7_thompson[2,1] = temp_Omega_hat_7_thompson[2,1] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) * (matrix(      1*U_hat[,t],N,1)) )
  # temp_Omega_hat_7_thompson[1,2] = temp_Omega_hat_7_thompson[1,2] + sum( matrix(      1*U_hat[,t-j],N,1) * (matrix(  X[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,2] = temp_Omega_hat_7_thompson[2,2] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) * (matrix(  X[,t]*U_hat[,t],N,1)) )
 }
 temp_Omega_hat_7_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_7_thompson
 Omega_hat_7_thompson = Omega_hat_7_thompson + temp_Omega_hat_7_thompson
}
Omega_hat_7_thompson = min(c(N,T))/T * Omega_hat_7_thompson

### COMPUTE M HAT ######################################################
S = colSums(X*U_hat)[2:T]
Slag = colSums(X*U_hat)[1:(T-1)]
rho_hat = cov(S,Slag)/var(Slag)
m = 1.8171 * (rho_hat^2/(1-rho_hat^2)^2)^(1/3) * T^(1/3)
m = (min(c(m,T-1)))

### OMEGA HAT 4TH TERM WITH BARTLET KERNEL #############################
Omega_hat_4 = matrix(0,2,2)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_4 = matrix(0,2,2)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_4[1,1] = 0 # temp_Omega_hat_4[1,1] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,1] = 0 # temp_Omega_hat_4[2,1] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[1,2] = 0 # temp_Omega_hat_4[1,2] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,2] = temp_Omega_hat_4[2,2] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
  }
  temp_Omega_hat_4 = w / (N^2*T) * temp_Omega_hat_4
  Omega_hat_4 = Omega_hat_4 + temp_Omega_hat_4
 }
 Omega_hat_4 = min(c(N,T))/T * Omega_hat_4
}

### OMEGA HAT 5TH TERM WITH BARTLET KERNEL #############################
Omega_hat_5 = matrix(0,2,2)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_5 = matrix(0,2,2)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_5[1,1] = 0 # temp_Omega_hat_5[1,1] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,1] = 0 # temp_Omega_hat_5[2,1] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
   temp_Omega_hat_5[1,2] = 0 # temp_Omega_hat_5[1,2] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,2] = temp_Omega_hat_5[2,2] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
  }
  temp_Omega_hat_5 = w / (N^2*T) * temp_Omega_hat_5
  Omega_hat_5 = Omega_hat_5 + temp_Omega_hat_5
 }
 Omega_hat_5 = min(c(N,T))/T * Omega_hat_5
}

########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################



########################################################################
# VARIANCE ESTIMATION (BEGIN)
########################################################################

### DESIGN MATRIX Q HAT ################################################
Q_hat = t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) / (N*T)

### STANDARD ERRORS ####################################################
SE_IID = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_iid )                                                                                                                       )%*%solve(Q_hat)/(N*T) ) )
SE_CRi = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CRt = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_2 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CGM = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 )                                                                                             )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_DDG = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 )                                                                                                           )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_THO = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4_thompson + Omega_hat_5_thompson - Omega_hat_6_thompson - Omega_hat_7_thompson ) )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CHS = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4 + Omega_hat_5 )                                                                 )%*%solve(Q_hat)/min(c(N,T)) ) )

########################################################################
# VARIANCE ESTIMATION (ENDED)
########################################################################



########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (BEGIN)
########################################################################
MNW_t_list = NULL
for( boot_iterations in 1:num_boot_MNW ){
rademacher = matrix(-1+2*(runif(N*T)<0.5),N,T)
U_hat_wild = rademacher * U_hat
Y_boot = beta_hat[1,1] + beta_hat[2,1] * X + U_hat_wild
beta_hat_boot = solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% t(cbind(1,matrix(c(X),,1))) %*% matrix(c(Y_boot),,1)
U_hat_boot = Y_boot - beta_hat_boot[1,1] - beta_hat_boot[2,1] * X

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1_boot = matrix(0,2,2)
for( i in 1:N ){
 # Omega_hat_1_boot[1,1] = Omega_hat_1_boot[1,1] + sum( t(matrix(    1*U_hat_boot[i,],1,T)) %*% matrix(    1*U_hat_boot[i,],1,T) )
 # Omega_hat_1_boot[2,1] = Omega_hat_1_boot[2,1] + sum( t(matrix(X[i,]*U_hat_boot[i,],1,T)) %*% matrix(    1*U_hat_boot[i,],1,T) )
 # Omega_hat_1_boot[1,2] = Omega_hat_1_boot[1,2] + sum( t(matrix(    1*U_hat_boot[i,],1,T)) %*% matrix(X[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[2,2] = Omega_hat_1_boot[2,2] + sum( t(matrix(X[i,]*U_hat_boot[i,],1,T)) %*% matrix(X[i,]*U_hat_boot[i,],1,T) )
}
Omega_hat_1_boot = min(c(N,T))/N/(N*T^2) * Omega_hat_1_boot

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2_boot = matrix(0,2,2)
for( t in 1:T ){
 # Omega_hat_2_boot[1,1] = Omega_hat_2_boot[1,1] + sum( matrix(    1*U_hat_boot[,t],N,1) %*% t(matrix(    1*U_hat_boot[,t],N,1)) )
 # Omega_hat_2_boot[2,1] = Omega_hat_2_boot[2,1] + sum( matrix(X[,t]*U_hat_boot[,t],N,1) %*% t(matrix(    1*U_hat_boot[,t],N,1)) )
 # Omega_hat_2_boot[1,2] = Omega_hat_2_boot[1,2] + sum( matrix(    1*U_hat_boot[,t],N,1) %*% t(matrix(X[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[2,2] = Omega_hat_2_boot[2,2] + sum( matrix(X[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X[,t]*U_hat_boot[,t],N,1)) )
}
Omega_hat_2_boot = min(c(N,T))/T/(N^2*T) * Omega_hat_2_boot

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3_boot = matrix(0,2,2)
# Omega_hat_3_boot[1,1] = t(matrix(c(1*U_hat_boot),,1)) %*% matrix(1*c(1*U_hat_boot),,1)
# Omega_hat_3_boot[2,1] = t(matrix(c(X*U_hat_boot),,1)) %*% matrix(1*c(1*U_hat_boot),,1)
# Omega_hat_3_boot[1,2] = t(matrix(c(1*U_hat_boot),,1)) %*% matrix(1*c(X*U_hat_boot),,1)
Omega_hat_3_boot[2,2] = t(matrix(c(X*U_hat_boot),,1)) %*% matrix(1*c(X*U_hat_boot),,1)
Omega_hat_3_boot = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3_boot

SE_CGM_boot = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1_boot + Omega_hat_2_boot - Omega_hat_3_boot ) )%*%solve(Q_hat)/min(c(N,T)) ) )
t_wild_boot = (beta_hat_boot - beta_hat ) / matrix(SE_CGM_boot,2,1)
MNW_t_list = cbind( MNW_t_list, t_wild_boot )
}
crit_MNW = apply(abs(MNW_t_list),1,quantile,p=0.95)
########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (ENDED)
########################################################################



cov_IID_95 = t(abs(beta_hat - true_beta) <= SE_IID*qnorm(0.975)); if( is.na(cov_IID_95)[1] ){cov_IID_95[1]=0}; if( is.na(cov_IID_95)[2] ){cov_IID_95[2]=0};
cov_CRi_95 = t(abs(beta_hat - true_beta) <= SE_CRi*qnorm(0.975)); if( is.na(cov_CRi_95)[1] ){cov_CRi_95[1]=0}; if( is.na(cov_CRi_95)[2] ){cov_CRi_95[2]=0};
cov_CRt_95 = t(abs(beta_hat - true_beta) <= SE_CRt*qnorm(0.975)); if( is.na(cov_CRt_95)[1] ){cov_CRt_95[1]=0}; if( is.na(cov_CRt_95)[2] ){cov_CRt_95[2]=0};
cov_CGM_95 = t(abs(beta_hat - true_beta) <= SE_CGM*qnorm(0.975)); if( is.na(cov_CGM_95)[1] ){cov_CGM_95[1]=0}; if( is.na(cov_CGM_95)[2] ){cov_CGM_95[2]=0};
cov_MNW_95 = t(abs(beta_hat - true_beta) <= SE_CGM*crit_MNW);     if( is.na(cov_CGM_95)[1] ){cov_CGM_95[1]=0}; if( is.na(cov_CGM_95)[2] ){cov_CGM_95[2]=0};
cov_THO_95 = t(abs(beta_hat - true_beta) <= SE_THO*qnorm(0.975)); if( is.na(cov_THO_95)[1] ){cov_THO_95[1]=0}; if( is.na(cov_THO_95)[2] ){cov_THO_95[2]=0};
cov_DDG_95 = t(abs(beta_hat - true_beta) <= SE_DDG*qnorm(0.975)); if( is.na(cov_DDG_95)[1] ){cov_DDG_95[1]=0}; if( is.na(cov_DDG_95)[2] ){cov_DDG_95[2]=0};
cov_CHS_95 = t(abs(beta_hat - true_beta) <= SE_CHS*qnorm(0.975)); if( is.na(cov_CHS_95)[1] ){cov_CHS_95[1]=0}; if( is.na(cov_CHS_95)[2] ){cov_CHS_95[2]=0};

list_IID_95 = rbind(list_IID_95, cov_IID_95)
list_CRi_95 = rbind(list_CRi_95, cov_CRi_95)
list_CRt_95 = rbind(list_CRt_95, cov_CRt_95)
list_CGM_95 = rbind(list_CGM_95, cov_CGM_95)
list_MNW_95 = rbind(list_MNW_95, cov_MNW_95)
list_THO_95 = rbind(list_THO_95, cov_THO_95)
list_DDG_95 = rbind(list_DDG_95, cov_DDG_95)
list_CHS_95 = rbind(list_CHS_95, cov_CHS_95)

options(digits=3)
print(c(iter,colMeans(list_IID_95)[2],colMeans(list_CRi_95)[2],colMeans(list_CRt_95)[2],colMeans(list_CGM_95)[2],colMeans(list_MNW_95)[2],colMeans(list_MEN_95)[2],colMeans(list_ME2_95)[2],colMeans(list_THO_95)[2],colMeans(list_CHS_95)[2]))
Sys.sleep(0.00000000000000001); flush.console()
}
########################################################################
# END OF MONTE CARLO ###################################################
########################################################################


RESULTS = cbind(colMeans(list_IID_95)[2],colMeans(list_CRi_95)[2],colMeans(list_CRt_95)[2],colMeans(list_CGM_95)[2],colMeans(list_MNW_95)[2],colMeans(list_MEN_95)[2],colMeans(list_ME2_95)[2],colMeans(list_THO_95)[2],colMeans(list_CHS_95)[2])
colnames(RESULTS) <- c("EHW_95","CRi_95","CRt_95","CGM_95","MNW_95","Men_95","Me2_95","Tho_95","CHS_95")
rownames(RESULTS) <- c("beta1")
RESULTS
if(DGP_AR){
 filename = paste("results_AR_f_",b[1],"_",b[2],"_",b[3],"_",b[4],"_",b[5],"_rho",rho,"_N",N,"_T",T,".csv",sep="")
}else{
 filename = paste("results_MA_f_",b[1],"_",b[2],"_",b[3],"_",b[4],"_",b[5],"_q",q,"_N",N,"_T",T,".csv",sep="")
}
write.csv(RESULTS,filename)
